#include "gc.h"
#include "itkImageRegionConstIteratorWithIndex.h"

void gen_area_file (const char* labelImagePathTemplate, 
		    const char* areaFilePathTemplate, 
		    int sliceNo0, int sliceNo1)
{
  char path0[200];
  char path1[200];
  sprintf(path0, labelImagePathTemplate, sliceNo0);
  sprintf(path1, labelImagePathTemplate, sliceNo1);
  LabelImage::Pointer labelImage0 = readLabelImage(path0);
  LabelImage::Pointer labelImage1 = readLabelImage(path1);
  LabelPixel maxLabel0 = 0;
  LabelPixel maxLabel1 = 0;
  itk::ImageRegionConstIteratorWithIndex<LabelImage> 
    itr0(labelImage0, labelImage0->GetLargestPossibleRegion());
  while (!itr0.IsAtEnd()) {
    if (itr0.Get() > maxLabel0) {
      maxLabel0 = itr0.Get();
    }
    ++itr0;
  }
  itk::ImageRegionConstIteratorWithIndex<LabelImage> 
    itr1(labelImage1, labelImage1->GetLargestPossibleRegion());
  while (!itr1.IsAtEnd()) {
    if (itr1.Get() > maxLabel1) {
      maxLabel1 = itr1.Get();
    }
    ++itr1;
  }
  int componentNum0 = maxLabel0 + 1;
  int componentNum1 = maxLabel1 + 1;
  int totalLength = componentNum0 * componentNum1;
  int* areaMatrix = new int[totalLength];
  for (int cnt = 0; cnt < totalLength; cnt++) {
    areaMatrix[cnt] = 0;
  }
  itk::ImageRegionConstIteratorWithIndex<LabelImage> 
    itr(labelImage0, labelImage0->GetLargestPossibleRegion());
  while (!itr.IsAtEnd()) {
    LabelImage::IndexType index = itr.GetIndex();
    LabelPixel label0 = itr.Get();
    LabelPixel label1 = labelImage1->GetPixel(index);
    areaMatrix[label0 * componentNum1 + label1]++;
    ++itr;
  }
  char path[200];
  sprintf(path, areaFilePathTemplate, sliceNo0, sliceNo1);
  writeArrayBinary(areaMatrix, totalLength, path);
  delete[] areaMatrix;
}



int main (int argc, char* argv[])
{
  clock_t startTime = clock();
  int argNo = 1;
  const char* labelImagePathTemplate = argv[argNo++];
  const char* areaFilePathTemplate = argv[argNo++];
  int startSliceNo = atoi(argv[argNo++]);
  int endSliceNo = atoi(argv[argNo++]);
  int sliceNoInc = atoi(argv[argNo++]);
#pragma omp parrallel for
  for (int sliceNo = startSliceNo; sliceNo <= endSliceNo - sliceNoInc; sliceNo++) {
    gen_area_file(labelImagePathTemplate, areaFilePathTemplate, 
  		  sliceNo, sliceNo + sliceNoInc);
    std::cout << "Slice " << sliceNo << " and " 
  	      << sliceNo + sliceNoInc << ": done." << std::endl;
  }
  std::cout << "Total time cost: "
	    << (double)(clock() - startTime) / CLOCKS_PER_SEC 
	    << " sec." << std::endl;
  return 0;
}
